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Abstract 

Electron transfer processes in Debye solvents are studied using a spectral 
analysis method recently proposed. Spectral structure of a nonadiabatic two- 
state diffusion equation is investigated to reveal various kinetic regimes char- 
acterized by a broad range of physical parameters; electronic coupling, energy 
bias, reorganization energy, and solvent relaxation rate. Within this unified 
framework, several kinetic behaviors of the electron transfer kinetics, includ- 
ing adiabatic Rabi oscillation, crossover from the nonadiabatic to adiabatic 
limits, transition from the incoherent to coherent kinetic limits, and dynamic 
bath effect, are demonstrated and compared with results from previous the- 
oretical models. Dynamics of the electron transfer system is also calculated 
with the spectral analysis method. It is pointed out that in the large reorgani- 
zation energy case the nonadiabatic diffusion equation exhibits a non-physical 
behavior, yielding a negative eigenvalue. 
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I. INTRODUCTION 

Since the seminal work of Marcus on the nonadiabatic electron transfer reactionsJlH a 
great amount of effort has been made in the studies of the electron transfer reaction! with 
a variety of tools such as time-resolved spectroscopies,IH§ computer simulation methodsjlh! 
and analytical theories.SHU It is not surprising that there have been great interests in the 
studies of the electron transfer reaction, considering that it is involved in many important 



chemical and biological systems (for the most recent reviews, see Ref. |23|). For example, 
in the photo-synthetic reaction center, electron transfer process creates the initial charge 
separation which will eventually lead to the production of the adenosine triphosphate.il 
Also, recent studies of the molecular electronics will depend crucially on the complete un- 
derstanding and controlling of the electron transfer in the chemical systems.il'il Another 
kind of electron transfer reactions which are currently subject to extensive studies is proton- 
coupled electron transfer reactions, and many theoretical@Hi3 and experimentalllHIil studies 
have been performed on this subject. 

One of important and ubiquitous aspects of electron transfer kinetics is the dynamic 
solvent effect on the electron transfer rate.&~3 As both experimental and theoretical in- 
vestigations have been carried out on electron transfer reactions in solutions, many diverse 
phenomena including predictions of the original Marcus theory have been revealed, depend- 
ing on physical parameters involved in the electron transfer kinetics.0~0&~S 

Electronic coupling constant, V, given by the interaction matrix element between the 
electron donor and acceptor wavefunctions is one of the most important physical parameters, 
and depending on its magnitude compared with other parameters, various kinetic regimes 
are exhibited in the electron transfer process. When the electronic coupling constant is 
the smallest parameter of the electron transfer process, the electron transfer rate is well 
described by perturbation theory, which predicts the golden-rule ratej! 

k ~ &gr = —ft—Pc (1-1) 
where p c is the equilibrium population of the reactant state in the crossing regime, and this 



is a well known result of the Marcus electron transfer theory. When the electronic coupling 
constant large enough, the overall reaction process is not determined by the Marcus rate, 
Eq. fll.ip , but by the solvent diffusion rate describing the polarization dynamics of the solvent 
molecules,^ 

k^k D ^ n\p c , (1.2) 

where Q is the solvent relaxation rate and A is the classical reorganization energy, and this 
case is called the solvent- controlled limit^® 

About two decades ago, Zusman0 investigated the crossover between the Marcus and 
solvent-controlled regimes in the studies of the electron transfer reaction in Debye solvents. 
He used a mixed quantum-classical approach where the stochastic operator is introduced 
to describe bath relaxation processes occurring on two diabatic surfacesffl§ and treated 
the electron transfer process as a nonadiabatic transition between them. Zusman solved the 
nonadiabatic diffusion equation in the weak coupling limit (V <C k B T), and obtained the 
expression for the overall electron transfer rate, 

k 1 = &GR + > (1-3) 

which shows a transition from the Marcus to the solvent-controlled limits in the nonadiabatic 
regime. 

As the electronic coupling constant is increased further to be comparable to or larger than 
the thermal energy, it is expected that the electron transfer process involves an adiabatic 
barrier crossing. Finally, when the electronic coupling constant is even larger so that it 
has the same order of magnitude as the solvent reorganization energy, A, which is indeed 
the case for the mixed valence compounds,0~0 the electronic states are delocalized on the 
lower adiabatic surface. Due to the derealization nature of electronic states, an adiabatic 
picture is more useful than the diabatic one for analyzing the short-time dynamics in strongly 
coupled systems.^! In this picture, electronic coherence arises from Rabi oscillation between 
two adiabatic surfaces. 



Although there have been several studies to bridge between the Marcus regime and 
the solvent-controlled regime using various approachesJlHl!'il~§l few studies have discussed 
diverse kinetic regimes in a unified way, and often different approaches are taken in different 
regimes. It is desirable to investigate the effects of solvent dynamics on the electron transfer 
process in a unified approach for various parameter regimes. 

As a general approach to describing condensed phase dynamics, we recently proposed a 
spectral analysis method^^ Instead of focusing on dynamical trajectories of the reduced 
density matrix for dissipative systems, this methodology investigates the spectral structure of 
the evolution operator for dissipative systems, and it has been applied to the electron transfer 
process in mixed valence compounds to investigate the possibility of electronic coherence in 
those systems.El 

In this paper we present a thorough analysis of the electron transfer kinetics in Debye 
solvents based on the spectral analysis method. Electron transfer rate constant extracted 
from the spectral analysis is compared with other previous results both in the nonadiabatic 
and adiabatic regimes, and transition from the incoherent to coherent regimes is demon- 
strated by the spectral analysis method. When the solvent relaxation rate is very fast, it is 
found that the solvent dynamics has a significant effect on the Marcus curve. The spectral 
analysis method is also utilized as a density matrix propagation scheme. Preliminary results 
of the spectral analysis method focusing on symmetric reaction cases have been reportedil. 

The rest of the paper is organized as follows. The spectral analysis method of the 
nonadiabatic diffusion equation is formulated in Sec. [HJ. Two important limiting cases are 
discussed in Sec. |II[ In Sec. |V| comprehensive analysis of the spectral structure of the 
nonadiabatic diffusion equation is performed for a broad range of parameters, and diverse 
kinetic behaviors in electron transfer reactions are identified and characterized. We conclude 
in Sec. |V| by summarizing results obtained in this work. 
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II. THEORY 



We consider two electronic states, |1) and |2), which represent the electron donor and 
acceptor sites of the electron transfer system, respectively, and they are coupled to each other 
via the electronic coupling matrix element, V. Moreover, each electronic state is coupled to 
bath degrees of freedom. 

There have been extensive studies of the solvent effect on electron transfer dynamics in 
literature with various approaches.0~0lil~il Considering that electron transfer processes are 
usually probed at room temperature in polar solvents, one can treat the bath degrees of free- 
dom classically. Zusman and Yakobson- Burshteinffl proposed a mixed quantum-classical 
evolution equation of the reduced density matrix, p(E,t), independently, to investigate the 
solvent effect on electron transfer, 

^p(E, t) = Cp(E, t) = (C B + iV)p(E, t). (2.1) 

Here, E is the solvent polarization energy which plays a role of the reaction coordinate as 
first noticed by Marcus.0 £ and V represent operators for the solvent relaxation dynamics 
and for the electronic transition between two states, respectively. Explicitly, Eq. (|2.1|) is 
written in terms of the density matrix elements, 

p n = Cnpn +iV(p 12 -p 2 i), (2-2) 

P22 = A2P22 ~ iV (Pii ~ P21), (2.3) 

V 

P12 = A2P12 - iu} 12 pi2 + i-r{pu - P22), (2.4) 

n 

V 

P21 = C21P21 + iv-upz! - z ^-(Pn - P22), (2.5) 

Here, the diagonal and off-diagonal matrix elements of the reduced density matrix p(E, t) 
represent populations of the electronic states and coherences between them, respectively, 
and Cij's describe the relaxation process of classical bath over the free energy surfaces, with 
Cu defined on the free energy surface for the zth electronic state, and £12 and £21 defined 
on the averaged free energy surface. The functional form for the free energy surface in the 
electron transfer system is usually harmonic, 



UiW = (2.6) 
U 2 (E) = + e : (2.7) 

where A is the reorganization energy, and we assume e < without loss of generality. It is 
convenient to define U and ThiO\ 2 are the average and the difference of the two free energy 
surfaces, respectively, 

U 1 (E) + U 2 (E ) = E 1 + * L 
K ' 2 4A 2' K J 

huj 12 {E) = U l {E)-U 2 {E) = E-e. (2.9) 

We set h — 1 for simplicity henceforth. This set of a mixed quantum-classical two-state 
equation has been previously derived by several authors§@§ starting from the spin-boson 
Hamiltonian. 

We note that many chemically and biologically important electron transfer processes 
take place in an overdamped solvent environment. Then, the bath relaxation operators in 
Eqs. ( |2.2] )-( PT5D are modeled by one-dimensional Fokker-Planck operators £y, 

where (3 = l/(k B T). 

The Fokker-Planck equation models the relaxation process of the density matrix element 
as a diffusion process in the energy space and various parameters are identified; energy 
diffusion constant, D-^ = QA 2 , fluctuation of the solvent polarization energy, A 2 = (E 2 ) = 
2\kftT, and characteristic timescale of a Debye solvent = I/O. The correlation function 
of the solvent polarization energy is given by a single exponential form in a Debye solvent, 

C(t) = (E(t)E(0)) = A 2 exp(-Ot). (2.12) 

Note that since the nuclear dynamics is modeled by the Fokker-Planck operator, the pos- 
sibility of the vibrational coherence is not considered in this model of electron transfer 
dynamics. 
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We investigate the spectral structure of the nonadiabatic diffusion operator, £, by cal- 
culating eigenvalues, {Z u }, and corresponding right and left eigenf unctions, an d 

C\^) = -Z v \i/£), (2.13) 
(tf\jC=-Z v (tf\. (2-14) 

Because the nonadiabatic Liouville operator is non-Hermitian, the eigenvalues are generally 

given by complex values, and the right and left eigenfunctions corresponding to the same 

eigenvalue are not simply the Hermitian conjugate to each other .0 

The method of eigenfunction solution is well known for the diffusion process on the 

i 



harmonic potential energy surface as discussed in App. Unlike the diffusion problem 
on the single potential energy surface, however, there have been limited studies on the 
nonadiabatic diffusion problem involving more than a single potential energy surface. Cukier 
and co-workers have calculated the electron transfer rate by calculating the lowest eigenvalue 
of the nonadiabatic diffusion equation in the weak-coupling regime.0 

In this paper, the eigenfunctions of C\2 are used as our basis set to represent the nona- 
diabatic diffusion equation. In principle, one could have chosen the eigenfunctions of C\\ or 
£22 as basis functions; however, in that case one has to evaluate appropriate Franck-Condon 
factors when calculating the coupling matrix elements. The Fokker-Planck operator £ 12 is 
defined on the averaged single harmonic potential centered at E — 0, and its eigenvalue 
solutions are obtained following a similar procedure in App. 

C l2 \(j>n) = -nn|0?>, (2.15) 
(rt\£ 12 = -n{l(rt\, (2.16) 

where n = 0, 1, 2, • • •, and the nth right and left eigenfunctions for C12 are given by 

<fa(E) = ^ -H n ( -ji- ) , (2.18) 

nV ; (2"n!)i(27rA 2 )3 \V2A J V ' 
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with H n being the nth order Hermite polynomial. 

We separate real and imaginary parts of the coherence density matrix, namely, u =Repi2 



and v =Impi2, and rewrite Eqs. (|2 . 2|) - (|2~5|) as 



pn = {C 12 + 5C)p n -2Vv, (2.19) 

P22 = {C12 ~ 5C)p 2 2 + 2Vv, (2.20) 

u = C 12 u + u 12 v, (2.21) 

v = C 12 v -u 12 u + V(p u -P22), (2-22) 



where 5C = (Cn — C 22 )/2. Then, all the relevant operators in Eqs. (|2.19|) - (|2.22|) can be 
evaluated in terms of the basis functions, 



(CIA2IC) = -nn8 nm , (2.23) 
(#|&C|0£> = _n^^V^ + T<W> (2.24) 



(</>£|u; 12 K) = ^2\k B T(yMSn,m~i + Vm + 16 n , m +i) - e^ nm , (2.25) 
(<t>n\V\<f>«) = V5 nm . (2.26) 

With this basis set, we can expand the density matrix elements as 

00 

Pll (E,t) = J2*n(t)<P*(E), (2.27) 

n=0 

00 

p 22 (E,t) = J2 b n(t)<f>«(E), (2.28) 

n=0 

00 

u(E 1 t) = y £c n {t)^(E), (2.29) 

n=0 

00 

v(E,t) = y Ed n {t)$(E!). (2.30) 



n=0 



Substituting Eqs. ( |2.27| )-( |2T30| ) into the right eigenvalue equation Eq. ( |2.13| ), we have the 



following coupled eigenvalue equations, 



- Z u a n = -nVta n - fiy ^ ^ \/na n _ 1 - 2Vd n , (2.31) 



-Z u b n = -nttbn + SI J -A-V^n-i + 2Vd n , (2.32) 
y 2k B T 



ZvC n = -nVLc n + J2\k B T(\/n + ld n+1 + ^Jnd n -i) - ed ni (2.33) 



Zydn = -nVtd n - J2Xk B T(Vn + lc n+ i + ^/nc n -i) + ec n + V(a n - b n ), (2.34) 



which is an explicit basis set representation for the two-state diffusion operator in Eqs. ( |2.2| )- 
Q2.5| ). Eigenvalue equations for the left eigenvector in Eq. ( |2.14| ) can be written by making 
the transpose of Eqs. ( |2.31| )-( |2.34| ). Diagonalizing the 4N x 4N matrix (N = number of basis 



functions) defined in Eqs. fl2.31|) -( p.34D , we obtain the eigenvalues Z u and the corresponding 



eigenvectors of the nonadiabatic diffusion operator, 

oo 

|^> = £iU0n>, (2-35) 

?1=0 

OO 

WI = E^I. (2-36) 

n=0 

where R nu and L vn are elements of the transformation matrices. 

In general, due to the non-Hermitian nature of the nonadiabatic diffusion operator, the 
left and right eigenfunctions do not form an orthogonal set by themselves. However, when 
the eigenvalues are all nondegenerate, the left and right eigenfunctions form an orthogonal 
and complete set in a dual Hilbert space.&^S Explicitly, we have 

oo 

L "nR nu > = 5 UU >, (2.37) 

n=0 
oo 

^ t RnvLyrn ^nrn: (2.38) 

for the orthogonality and the completeness relation, respectively. Using these properties, we 
can construct the real time propagator for the operator £ as 



G(t)=£lV£><V£|e- z "*, (2.39) 
and express the time evolution of the density matrix as an eigenfunction expansion, 

oo 

\p(t)) = G(t)\p(0)) = £ |^)(^|p(0))e-^. (2.40) 

Note that the right and left eigenfunctions play asymmetric roles in the construction of the 
propagator for the non-Hermitian operator. 



III. LIMITING CASES 



In order to compare the eigenvalue solution developed in this work with previous theo- 
retical predictions available in different kinetic regimes, we briefly discuss two limiting cases 
which have been studied extensively in the literature .0Hll Instead of giving detailed deriva- 
tions, we will briefly mention the solutions in these limiting cases, relegating details to the 
App. § 



A. Nonadiabatic Regime : Weak Coupling Case 



When the electronic coupling matrix element V is very small, we can reduce the full 
nonadiabatic diffusion equation, Eqs. (|2.2|)-(|2.5|), into the population evolution equation, 



and detailed derivations can be found in Refs. 16,17 (see also App. |B|). Kinetic equations 
for the population elements of p can be approximately written as 

dpn 



dt 
dt 



C u pn - K(p n - p 22 ), 
C22P22 - K(p 22 - P11), 



(3.1) 
(3.2) 



where K(E) is the rate kernel given in Eq. ( p5|) . A dynamical quantity usually measured 
in the electron transfer kinetics experiment is the total population, Pi(t), in each electronic 
state rather than the polarization energy dependent population, pa(E,t), 



Pi(t) 



dE Pi (E,t). 



(3.3) 



Using the projection operator method and making a time-scale separation approximation, 
it can be shown that the kinetic process between two electronic states is described by the 
time-independent rate constant instead of the rate kernel in the weak coupling limit,0'0( 
also App. O). Then the kinetic equation for Pi is given by00 



sec 



( p. i+\ \ ( p U p. i+\ \ 



Pl(t) 

A(t) 



k 1 -k 2 



Pi(t) 

P2(t) 



(3.4) 
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Here k 1 and k 2 are forward and backward rate constants and they are related to the nona- 
diabatic transition rates k^ A and solvent diffusion rates k^ 

k = - JT~[_ T~2 / ; 2 ' (3.5) 

1 + k NA /K D + Kna/^D 
^NA = ^NA ^NA' (3-6) 

&nd = k 1 + k 2 = k J A , (3.7) 

1 -+- %A/ K D ~+~ %A/ K T> 

Recently, Eq. (|3.7|) has also been obtained for the symmetric reaction case in the weak cou- 
pling limit by calculating the first excited eigenvalue explicitly via the Goldstone theorem.il 
Eq. ( |3.7| ) has a form of an overall relaxation rate for consecutive reactions, and it involves 
two different types of rate processes; nonadiabatic transition and solvent diffusion rate. 
Nonadiabatic transition rate constants k^ A and k^ A2 are forward and backward quantum 
transition rates between the two electronic states and are calculated in terms of the coherent 
Green's functions, G 12 {E,t\E ) and G 2 i{E, t|£ )ffl(see also App. fj), 

roc 

fc* = 27rV 2 Re / dte l(e+x)t - 9{t \ (3.8) 
Jo 

where the bath correlation function g(t) is given by 



g(t) 



[exp(-fit) + ttt - 1] . (3.9) 



with +(— ) sign for i = 1(2). When the bath dynamics is slow such that Q -C A, Eq. 
reduces to the standard Marcus result (see App. |B|), 

^A-^R^TTVVrOO, (3.10) 

fc NA ~^ R + 4 R . (3.11) 



where Pi q (e) = e~^ ±£ ) / 2A / y/27rA 2 is the equilibrium population of the ith state at the 
crossing point E = e. When Q ~ A, the nonadiabatic transition rate shows a deviation 
from the Marcus rate, which will be discussed later in details. 

Solvent diffusion rate constants, and k 2 ^, describe solvent relaxation processes that 
equilibrate the nonequilibrium wavepacket created at the crossing region, E c = e, during the 
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electron transfer process. They can be written in terms of the population Green's functions 
G n (E,t\E ) and G 22 (E, t|£ o ),0ffl 

^~ 1 = fHejJo d ^Me)-G«(e,oo|e)]. ( 3 - 12 ) 

Noting that Gu(e, oo|e) = Pi q (e) is given by the equilibrium population distribution of the 
zth state at E = e, the solvent diffusion rates are identified as the inverse of the mean survival 
time of the relative nonequilibrium population created at the crossing point, E c = e. 

The boundary between the Marcus regime and the solvent-controlled regime in the nona- 
diabatic limit is estimated when comparing the Marcus rate in Eq. ( |3.10| ), and the solvent- 
diffusion rate in Eq. ( |B13| ), which leads to 

2ttV 2 ~ |e± A|a (3.13) 

The nonadiabatic regime is established when the timescale of the off-diagonal density matrix 
is much faster than that of the diagonal ones, yielding the following condition,0H 

V < D^ 3 = {2Xk B TQ) 1/3 . (3.14) 

in addition to the weak coupling condition, V <C k B T, for the validity of the nonadiabatic 
regime. 



B. Adiabatic Regime : Strong Coupling Case 



When the electronic coupling constant is comparable to or larger than the thermal energy 
j3V ^ 1, it is more adequate to describe electron transfer in an adiabatic representation. 
Given a specific bath configuration characterized by the value of E, the electronic part of 
the Hamiltonian in the diabatic representation is given by 



H A (E) 



1 ' U^E) V ^ 



(3.15) 



/ 



V U 2 (E) 

After the diagonalization of H&(E) for a given bath configuration E, we can obtain two 
adiabatic surfaces, 
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U±{E) = U(E) ±^[uj 12 (E) 2 + 4V 2 f /2 . (3.16) 

The separation between two adiabatic surfaces is much larger than ksT, U + — Z7_ > 
V ^> ksT. When the energy bias e and the electronic coupling V are not so large compared 
with the reorganization energy, A, the lower adiabatic surface, £/_, has a well-defined double 
well structure. In this case, we can consider the electron transfer process as a diffusional 
barrier crossing process occurring on the lower adiabatic surface Dad = U- described by the 
following diffusion equation in the energy space, 

% = *{m +fi m v 'A'' =c »»- (3 - 17) 



To calculate the reaction rate for the adiabatic barrier crossing described by Eq. ( |3.17| ) 



one may take several different routes. For example, one can calculate the barrier crossing 
rate as the first nonzero eigenvalue of the Fokker-Planck operator, £ad> corresponding to 
the lower adiabatic surface,^ 

£adV> = -Zip fc A v D = Z x (3.18) 

Another popular approach is a population-over-flux method,i3 which obtains the barrier 
crossing rate by calculating the steady-state flux at the crossing region, E c = e, divided by 
equilibrium population in the reactant region, 



[ Ec dE 1 r dE 2 e^ u -^- u ^ 

J—oo J —oo 



-1 



(3.19) 



In the symmetric reaction case (e = 0), the adiabatic rate given in Eq. ( |3.19| ) can be evaluated 



approximately in the strong coupling limit (f3V ^> 1) by expanding the integrand at E = 
and performing a Gaussian integration approximately to yield0'@ 



k ™ K i;m e ~^~ v) ' (3 - 20) 



which is the Kramer's adiabatic reaction rate in the strong damping regime.0 
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IV. RESULTS OF SPECTRAL ANALYSIS 



A. General Feature of Spectra 

Figure [I] shows the eigenvalue spectra in three different cases of the energy bias for 
various electronic couplings. Since the nonadiabatic diffusion operator is not Hermitian, the 
eigenvalues are complex in general, 

&l>„ = -Z^ u = -{Z' v + %Z'l)^ v . (4.1) 

In order for the density matrix relax to the stationary state at long times, the real part of 
the eigenvalue should satisfy Z' v > 0. Eigenvalues of the nonadiabatic diffusion operator 
are classified into three different cases with their associated kinetic behaviors in the density 
matrix evolution; 

(1) Z' Q = Zq = : equilibrium, 

(2) Z' v > 0, Z" = : exponential decay, 

(3) Z' v > 0, Z" : damped oscillation. 

The eigenstate with zero eigenvalue [case (1)] corresponds to the equilibrium solution of the 
density matrix while those with real [case (2)] and complex eigenvalues [case (3)] correspond 
to exponential decays and to damped oscillations in the evolution of the density matrix, 
respectively. When the first nonzero real eigenvalue Z[ is well separated from the other 
higher eigenvalues, Z[ <C Z' u>2 , the dynamics of the density matrix will be described by a 
single exponential process, and the relaxation rate is well defined as k = Z[. However, if this 
is not the case, the dynamics of the density matrix will generally involve multiple timescales 
and coherent oscillations. 

We used 800 basis functions to calculate the eigenvalues, and to remove the finite basis set 
effect at higher eigenvalues, we present only lower 400 eigenvalues in Fig. [I]. Three different 
cases of the energy bias chosen in Fig. |l|(a)-(c) correspond to different regimes in nonadiabatic 
electron transfer theory, i.e., normal (|e| < A in Fig. 0(a) ), activationless (|e| = A in Fig. [j](b)), 
and inverted regimes (|e| > A in Fig. |l|(c)), respectively. As a general feature, all the 
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spectra show tree-like structures with three major characteristic branches; one branch for the 
real eigenvalues in the middle, real axis, and two other branches for the complex conjugate 
eigenvalues. Therefore, the spectral structure indicates that multiple exponential decays (real 
eigenvalues) and damped oscillations (complex eigenvalues) are inherent features in electron 
transfer processes in the overdamped solvent. In order to infer the dynamical behavior of 
the density matrix, we only need to focus on the lower part of the eigenvalue diagrams in 
Fig. 1 corresponding to the eigenstates with Z'/Q ^ 1 since those with Z'/Q ^> 1 will decay 
out very quickly. 

(i) symmetric reaction case(e = 0) : In Fig. 0(a), one can notice that the real parts of 
the eigenvalues, Z' u J s, located at the branching regime in the eigenvalue tree decrease as the 
coupling constant increases up to (3V < 1, which means that the oscillatory components 
will persist during a more extended period t ~ 1/Z'„ as t ne coupling constant increases. 
When the coupling constant increases further, (3V 3> 1, the branching point of the complex 
branches start to separates along the imaginary axis and the separation between the two 
complex branches along the imaginary axis is given by ~ 2V, and this corresponds to the 
Rabi oscillation frequency for the two adiabatic states in the strong coupling regime 

(ii) asymmetric reaction cases(e ^ 0) : In Figs. [TJ(b) (/?|e| = 10) and (c) (/3|e| = 20), we 
notice that spectral structures show more branched behaviors in the asymmetric reaction 
cases than in the symmetric reaction case in Fig. |l](a). When comparing the zero electronic 
coupling cases, (3V = 0, in Figs. |l|(a), (b), and (c), we notice that in the asymmetric 
reaction cases two complex branches separate out along the imaginary axis by the amount 
of the energy bias, e unlike in the symmetric reaction case. This corresponds to the natural 
oscillation frequency of the off-diagonal density matrix elements. 

In Fig. |2| we have followed the evolution of the lowest 30 eigenvalues in both the real and 
imaginary axes as a function of the coupling constant for different energy bias cases. In all the 
three cases real parts of the eigenvalues start as pairs of doubly degenerate states separated 
from each other by the solvent relaxation rate, /3Q = 1, when V = 0. As the coupling 
constant increases, the degeneracies of the real eigenvalues are first lifted at small coupling 
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constants, and bifurcation and coalescence of eigenvalues occur at large coupling constants. 
The coalescence corresponds to the case where two closely separated real eigenvalues become 
complex conjugate with the same Z', and this indicates the transition from the incoherent, 
rate process to coherent oscillatory behavior in the dynamics of the density matrix, which 
will be further discussed in Sec. |1V C\ In the imaginary part diagram, we find an interesting 
behavior that imaginary parts of the eigenvalues asymptotically follow the Rabi oscillation 
frequency, 

Z" « ^Rabi = vW 2 + e 2 , (4.2) 

in the strong coupling limit, and this demonstrates the signature of adiabatic Rabi oscillation 
picture of the electron transfer in the strong electronic coupling case. 



B. Crossover from Nonadiabatic to Adiabatic Regimes 

In Fig. |3] we compared the first nonzero eigenvalue of the nonadiabatic diffusion operator 
with several predictions from existing theories in order to characterize various kinetic regimes 
in the nonadiabatic diffusion equation. The filled circles are calculation results of the first 
nonzero eigenvalue obtained from the spectral analysis method. We compare our eigenvalue 
solution with previous theoretical results, each of which is applicable for different limits; 

(i) weak coupling limit : The eigenvalue solution is compared with the purely nonadi- 
abatic transition rate without the solvent diffusion effects, Eq. Q3.6Q (solid line) and with 
the nonadiabatic-diffusion rate with the solvent diffusion effects included, Eq. ( |3.7|) (long 
dashed line). In any case, we compare the eigenvalue solution with the sum of the forward 
and backward rates, since the eigenvalue solution gives the overall relaxation rate of the 
density matrix. 

(ii) strong coupling limit : We compare the eigenvalue solution with the adiabatic barrier 
crossing rates as discussed in Sec. [Ill B|. We calculated two different adiabatic reaction 



rates, the lowest nonzero eigenvalue given in Eq. ( ft.lSQ (dot-dashed line) and the flux-over- 
population rate given in Eq. (|3.19|) (short dashed line), respectively. 
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In the slow bath relaxation case (j3Q = 0.05) given in Fig. 0(a), when the the coupling 
constant is very small (J3V ^ 0.1) (Marcus regime), the eigenvalue solution agrees well 
both with the nonadiabatic transition rate and with the nonadiabatic-diffusion rate. The 
nonadiabatic transition is the rate limiting step in this case. As the coupling constant 
becomes large and so does the nonadiabatic transition rate, the overall rate is now affected 
by the solvent relaxation rate, and the eigenvalue solution follows the nonadiabatic-diffusion 
rate, Eq. (|3.7| ), and this is the solvent- controlled regime. 

As the coupling constant becomes much larger than the thermal energy, the eigenvalue 
solution shows a saturation behavior at f3V ~ 1 first, and starts to increase rapidly as the 
electronic coupling becomes stronger, and this indicates the adiabatic barrier crossing in the 
strong coupling regime. In this case, (3V ^> 1, the eigenvalue shows a qualitative agreement 
with the adiabatic reaction rates, Eq. (|3.18|) and Eq. (|3.19|) . 

As we increase the solvent relaxation rate (j3Q = 1) in Fig. 0(b), the distinction between 
the nonadiabatic and the adiabatic regime in the eigenvalue solution becomes less obvious 
than in the slow bath case, Fig. 0(a). In the fast bath relaxation the nonadiabatic 

transition rate becomes large to be in the solvent-controlled regime, the adiabatic barrier 
crossing is already appreciable, so the transition between the solvent-controlled and adiabatic 
barrier crossing regimes is not so clear. This is related to a dynamical modulation of the 



nonadiabatic transition rate in the fast bath case, and will be discussed in details in Sec. [IV D 
The agreement between the eigenvalue solution and the two adiabatic rate calculations 
is qualitative in the strong coupling case. Although it is reasonable to treat the electron 
transfer process as an adiabatic barrier crossing reaction on the lower adiabatic surface 
when the coupling constant is large (flV 3> 1), it should be mentioned that a rigorous 
mathematical proof has not been given to show that the electron transfer rate from the 
nonadiabatic diffusion equation really corresponds to the adiabatic reaction rate in the 
strong coupling limit. Therefore, it may not be surprising that the agreement between the 
eigenvalue solution and the adiabatic solutions is only qualitative. Related to this problem, 
it should be mentioned that the original nonadiabatic diffusion equation is based on diabatic 
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representations of electronic states, therefore it may not give quantitatively the same rate 
constant as that from the adiabatic diffusion equation in the strong coupling case. 

The transition from the Marcus to solvent-controlled regimes has been studied theoreti- 
cally by using different theoretical methods such as the projection operator techniques@'0ii 
and the path-integral methods.il These theoretical approaches do not cover the adiabatic 
barrier crossing regime. The crossover behavior from the nonadiabatic to adiabatic regimes 
has been observed in quantum mechanical approaches to electron transfer theories based on 
the instanton solution&~H and the diagrammatic technique^ which do not take into account 
the solvent dynamics effects. The eigenvalue solution presented here clearly demonstrates 
transitions between three different regimes. 

C. Transition from Incoherent to Coherent Regimes 

For the symmetric electron transfer reaction case(e = 0) the lower adiabatic surface, L7_, 
has a double well structure if V < A/2. When the electronic coupling constant becomes 
even larger than A/2, the lower adiabatic surface becomes a single well without any barrier. 
We consider this situation as a thermodynamic transition from the localized to delocalized 
electronic states when viewed from a perspective of the diabatic states. Even though the 
free energy surfaces of the donor and acceptor states support the large energy barrier be- 
tween them in the diabatic picture, a large coupling constant in the same order of magnitude 
as the reorganization energy makes the distinction between the donor and acceptor states 
inappropriate. Instead, two new single well potential surfaces are obtained, whose eigen- 
functions are neither donor nor acceptor wavef unctions, and rather linear combinations of 
them. This situation occurs in mixed valence compounds and other strongly coupled systems 
where V ~ A ^> /c^T.&B The critical value of V = V c at which the double well structure 
disappears in the lower adiabatic surface, U- can be obtained from Eq. fl3.16| ) when the 
reorganization energy is A and the energy bias is |e|, 

(2V C ) 2/3 + |e| 2/3 = A 2/3 , (4.3) 



which indicates that the localization-delocalization transition occurs at a smaller value 
of V c < A/2 in the asymmetric reaction case. Recently, Cukier and co-workers, and 
Hanggi and co-workers have studied such strong coupling regimes based on the semi-classical 
approachesJ^S 

To investigate the feature of the localization-delocalization transition in the spectral 
structure, we show the lowest 5 eigenvalues as a function of the coupling constant in the 
strong coupling regime when e = in Fig. |j. When the electronic coupling is about the half 
of the reorganization energy as predicted by Eq. (|4.3| ) f3V ^ 5, the first nonzero eigenvalue 
(solid line) is no longer well separated from higher complex eigenvalues, it is no longer valid 
to describe the dynamics of the density matrix by a incoherent rate process. Also, there is 
a rapid drop of the first excited eigenvalue around (5V ~ 6 due to the crossing of the first 
excited state eigenvalue and higher complex eigenvalues (note that (3V C = 5 in this case). 
The crossing behavior of the first and higher excited state eigenvalues is a signature of the 
transition between the incoherent to coherent regimes in electron transfer processes. 

A crossing of the two lowest nonzero eigenvalues can be used as an approximate criterion 
for the transition from the incoherent to coherent regimes. In Fig. (|^), we plotted the value 
of the electronic coupling constant, V c , where a crossing of the two lowest nonzero eigenvalues 
occur for a given energy bias, |e|, and the prediction of Eq. (|4.3| ) for the case (3VL = 0.01 
and j3X = 10. They show the same trend that the coherent regime appears earlier in the 
asymmetric reaction case than in the symmetric reaction case. 

D. Dynamic Bath Effect 

In addition to the fact that the overall electron transfer rate appears a combined result 
of the nonadiabatic transition and solvent diffusion processes in Eq. ( |3.7| ), the nonadiabatic 
transition rate itself, Eq. ( |3.8| ), is modulated by the dynamic nature of the solvent relaxation 
process that is characterized by the bath correlation function in Eq. ( |3.9| ). We will call this 
a dynamic bath effect. This will be evident if the bath relaxation rate is comparable to the 
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amount of the polarization energy fluctuation, 

n ~ A = ^2\k B T. (4.4) 

To investigate the dynamic bath effect, we compare the eigenvalue solution and two 
different reaction rates, Marcus rate in Eq. ( |3.f f| ) (static bath limit) and the nonadiabatic 
transition rate, Eq. (|3.6| ) (dynamic bath effect included) in the weak coupling limit as we vary 
the energy bias in Fig. || All the calculation results show a characteristic behavior of the 
Marcus curve; they show a maximum value when |e| ~ A. In the slow bath case, Fig. §(a), all 
the three curves agree with each other since the dynamic bath effect is negligible in this case. 
In the fast bath case, Fig. ||(b), we notice that the eigenvalue and the nonadiabatic transition 
rate agree with each other, and show a deviation from the Marcus rate. In particular, in 
the fast bath case, the maximum position of the eigenvalue as well as of the nonadiabatic 
transition rate shift toward a smaller value of the energy bias compared with the Marcus rate. 

In the static bath limit, the nonadiabatic transition rate is reduced to the Marcus rate, 
exhibiting the maximum position at the activationless case, |e| = A. When the bath re- 
laxation is fast, the solvent reorganization energy is effectively reduced, and the maximum 
electron transfer rate appears at a smaller value of the bare solvent reorganization energy, 
exhibiting the peak position in the normal regime, |e| < A, and this results in a significant 
deviation of the nonadiabatic transition rate from the Marcus rate. The amount of the peak 
shift as a function of the solvent relaxation rate can also be estimated from Eq. (|3.8| ), and 
it is shown in Fig. [7] that the peak shift increases linearly with Q. 

The peak shift observed in the eigenvalue solution has a close analogy to the Stokes 
shift observed in the condensed phase spectroscopy.!! It is worthwhile to mention that g(t) 
obtained in Eq. ( |3.9| ) has exactly the same form as the line broadening function obtained in 
the hight temperature limit of the overdamped Brownian oscillator model in the condensed 
phase spectroscopy]" and the nonadiabatic transition rates given in Eq. ( |3. 8| ) have the 
same functional forms as the line shapes of the chromophore interacting with the Debye 
solvent. As discussed in the App. ||, the Stokes shift between the absorption and the 
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emission profiles disappears in the fast bath case and it has the same origin as the peak shift 
observed here. 

Implication of the dynamic bath effect on the experimental observables may be signifi- 
cant. The reorganization energy of the solvent can be estimated from the maximum of the 
electron transfer rates based on a series of the measurements of electron transfer rates for 
different systems of which the energy bias and the electronic coupling constants are known. 
However, the dynamic bath effect can make the value of the reorganization energy measured 
in the experiment smaller than its true value. Therefore, care should be taken to disentangle 
the dynamical bath effect from the measured value of the solvent reorganization energy. 



E. Density Matrix Propagation 

The spectral method can be used as a propagation method of the density matrix as shown 
in Eq. ( [2.40| ). We demonstrate the time evolution of the donor and acceptor populations, 
Pi(t) and P 2 (t) based on the spectral method. Parameters are chosen as f3\ = 10, j3Q = 1, 
and f3e = 0. The electronic coupling constant will be varied. The initial condition is chosen 
such that at t = only the donor state is populated with the equilibrium distribution in the 
diabatic state, 

p u (E,0) = p%(E), (4.5) 
p 22 (E, 0) = p 12 (E, 0) = p 21 (E, 0) = 0. (4.6) 

This initial condition corresponds to that of the photo-induced back electron transfer 
react ion .EHEIS EI 

In Fig. |H we show the time evolution of the donor (solid line) and the acceptor (dashed 
line) state populations with various electronic coupling constants. At the smallest coupling 
constant, (3V = 0.5, the donor and the acceptor populations relax exponentially to the 
equilibrium at long times. The electronic coupling constant is small such that the time evo- 
lution is governed by the lowest nonzero eigenvalue of the nonadiabatic diffusion equation. 
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As the electronic coupling constant is increased, the temporal behavior of the density ma- 
trix undergoes the incoherent- coherent transition as discussed based on the spectral analysis 
in Sec. [IV D\ When the coherent oscillation is observed in the temporal evolution of the 
density matrix(/3V = 2 and (3V = 4), the oscillation frequency matches the Rabi oscillation 
frequency given in Eq. (|4.2|). The time-dependent study of the nonadiabatic diffusion equa- 
tion based on the spectral method was also performed for back electron transfer dynamics 
occurring in the mixed-valence system.0 



F. Nonphysical Behavior in Nonadiabatic Diffusion Equation 

We found that the nonadiabatic diffusion equation can yield nonphysical behavior in 
some parameter regimes. We checked the case of large reorganization energy, /?A> 1, and 
found that in this case the smallest eigenvalue can be negative instead of zero, which yields 
an exponential growth of the density matrix in the long time limit. Fig. ^] demonstrates 
how two lowest real eigenvalues change as the reorganization energy is varied in the weak 
coupling case. When f3\ > 12, the lowest nonzero eigenvalue becomes negative instead of 
positive, thus violating the positive defmiteness of the density matrix in the nonadiabatic 
diffusion equation. 

There are several possible reasons for this problematic behavior. One of them is an 
intrinsic limitation in consistent treatments of the electronic and nuclear degrees of freedom 
in the mixed quantum-classical approach. Similar problem in the mixed quantum-classical 
approach has been pointed out in other contextsSEll The other possible reason for the 
positivity violation may be the Markovian approximation used in the calculation of the 
frictional kernel in the overdamped regime for the nonadiabatic dynamics. Although the 
Markovian approximation for the friction kernel has been successfully applied to the adiabatic 
reaction rate theory where the single potential energy surface is involved,^ there have been 
limited studies on the validity of this approximation in the nonadiabatic reaction cases-SEi 
Along this line, similar nonphysical behaviors in the nonadiabatic diffusion equation yielding 
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negative electron transfer rates have been found recently-L^To The precise origin of this 
artifact is still not clear, and more studies need to be done to clarify this issue in a decisive 
way. 



V. CONCLUSION 

In this paper we studied the electron transfer reaction kinetics in Debye solvents described 
by the two-state nonadiabatic diffusion equation. By applying the spectral analysis method 
developed recently, we investigated various kinetic regimes of the electron transfer in Debye 
solvents in a unified way. The results obtained in this work are summarized; 

(i) The general spectral structure of the nonadiabatic diffusion equation has been re- 
vealed. It exhibits a tree-like structure with three major branches; a single branch of real 
eigenvalues corresponding to multi-exponential decays and two symmetric branches of com- 
plex conjugate eigenvalues corresponding to damped oscillations. 

(ii) Several kinetic regimes in the electron transfer reactions in solutions have been 
identified. The first nonzero eigenvalue solution bridges between the Marcus, the solvent- 
controlled, and the adiabatic crossing regime in the incoherent kinetics regime. It agrees 
well with nonadiabatic reaction rates from previous theoretical results in the weak cou- 
pling regime. In the adiabatic crossing regime the eigenvalue solution agrees with adiabatic 
reaction rates qualitatively. 

(iii) When the electronic coupling constant is about half the reorganization energy in the 
symmetric reaction case (less than half in the asymmetric reaction case), the eigenvalue dia- 
gram shows coalescence/bifurcation behavior, where the two lowest nonzero eigenvalues cross 
each other and become complex conjugate. This indicates the damped oscillation behavior 
in the overdamped solvent, and correlates with the localization-delocalization transition in 
the lower adiabatic surface. 

(iv) Due to the finite timescale of the solvation, the eigenvalue solution is modulated 
in the weak coupling regime such that the maximum position of the electron transfer rate 
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appears at a smaller value of the bare solvent reorganization energy in the Marcus curve, 
which shows a substantial deviation from the standard Marcus theory. The consequence of 
the dynamic bath effect may be significant in the experimental determination of the solvent 
reorganization energy. 

(v) The spectral analysis method was used as a propagation scheme to study the time- 
dependent behavior of the density matrix, and it has confirmed previous predictions made 
based on the investigation of the spectral structure. 

(vi) Nonphysical behavior of the nonadiabatic diffusion equation, violating the positivity, 
has been identified in the large reorganization energy case by the spectral analysis method, 
and the precise origin of this problem should deserve further investigations. 

In the isolated quantum system the eigenvalue solution method offers a powerful way 
to analyze the dynamics. In comparison, the spectral structure of the dissipative quantum 
dynamics has not been fully explored yet, and most studies on the quantum dissipative 
system employ trajectory methods. In this regard, our studies on the spectral structure 
of the electron transfer kinetics can be a good motivation to alternative studies of various 
dissipative systems. The nonadiabatic diffusion equation studied in this work is based on 
the mixed quantum-classical approach. However, the spectral analysis method itself can be 
applied to other quantum dissipative dynamical equations in general such as the Redfield 
equation,@'0 quantum Fokker-Planck equation,0E§ and quantum master equation,^! and the 
studies based upon the spectral analysis will bring more direct insights on the dissipative 
dynamics. 
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APPENDIX A: SPECTRAL ANALYSIS FOR FOKKER-PLANCK OPERATOR 
ON A SINGLE HARMONIC POTENTIAL 



For a single harmonic potential, U(x) = ^mu 2 x 2 , the Fokker-Planck operator is given by 

^ = D (& + B i u ')- < A1) 

and it is possible to transfrom Cpp into a quantum mechanical Hamiltonian corresponding 
to a fictitious harmonic oscillator in imaginary time, 

H = - e ^(*)/2£ Fpe -W*>/2 = -— ^- + V(x), (A2) 

2 fx ox 1 

where fx = (2D) -1 and the potential for the fictitious harmonic oscillator is given by 
V(x) = D 



\(PU'(x)f-l(3U"(x) 



\fi^ 2 x 2 - J, (A3) 



with 7 = Dmuj 2 /ksT . Since the transformed potential in Eq. ( |A3|) is in the same form as 
that of a harmonic oscillator model with an off-set of the zero point energy, eigenvalues and 
eigenfunctions for the original Fokker-Planck operator can be obtained immediately from 
the eigenvalue solutions of the harmonic oscillator Hamiltonian,ii 

C FP \<) = -njM, (A4) 
(V£|£fp = -nT^nl, (A5) 

where i/j%(x) = H n (x/ y/2x ) / (prfPynlxo) 1 / 2 and ifj*(x) = e~ x% l 2x ^{x) with x = 



IzbT jmuj 2 and H n being the nth order Hermite polynomial. 



APPENDIX B: RATE EQUATIONS IN NONADIABATIC REGIME 



We briefly review the rate expressions in the nonadiabatic limit {(3V -C l).ffl First, we 
obtain a formal solution for pi2(E,t) in the time domain from Eq. ( |2.4| ), 



pt poo 

p 12 (E, t)=iV dr dE G 12 (E, t\E ) \p n (E , t-r)- p 22 (E , t - r)] , (Bl) 

J J-oo 
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where Gu(E,t\E ) is the Green's function for the operator, £ i2 — iuJvz, 

G 12 (E,t\E ) = (E\e( c »-****\Eo). (B2) 



Substituting pi 2 and p 2 i = p\ 2 from Eq. (|B1|) into the equations for pn and p 22 in Eqs. ( [2~2l 
and (|2.3p , we can obtain a closed set of integro-differential equations for pn and p 22 , 



d 

—pu(E,t) = CiiPii{E,t) 

/oo /"t 
d£ / drGi 2 (E,r|S )[Pii(^o,t - r) - Pn(E ,t- r)}. (B3) 
-oo JO 

In the weak coupling regime, f3V -C 1, the electron transfer process is characterized as 
nonadiabatic, and the transition between two electronic states takes place only when the 
bath polarization energy is close to the free energy gap, and pn and p 22 vary slowly compared 
with pi2 and p 2 i. We can assume that the coherent Green's function G\ 2 varies much faster 
than the population difference pn — p 22 in this case and that dominant contribution to 
the integral in Eq. ( |B3| ) comes from the short time region. Using these arguments we can 
approximately deconvolute the integrals and extend the upper limit of the integral to the 
infinity, 

/oo rt 
dE / drG 12 (E,T\E )[p u (E ,t-r)-p 22 (E ,t-r)] 
-oo JO 

^K(E)[ Pu (E,t)-p 22 (E,t)}, (B4) 

where K(E) is the rate kernel and is given in terms of the coherent Green function G(E, t\Eo) 
defined in Eq. flB^JHB 

roo poo 

K(E) = 2V 2 Re / dr / dE G 12 {E,r\E ). (B5) 

JO J-oo 

The approximations made here are well established in the normal regime of the electron 
transfer, but they are questionable in the inverted or activationless regime, and improvements 
on these regimes have been made.0 

The nonadiabatic rates are given in terms of the coherent Green's function]^ 1 1 



roo roo roo 

k\, A = 2V 2 Re dt dE dE G tJ (E,t\E )pf(E ). (B6) 

JO J —oo J —oo 
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The coherent Green's function, Gi 2 (E,t\Eo), can be evaluated first by transforming the 
Fokker-Planck operator into the harmonic oscillator Hamiltonian as done in Eq. (|A2| ), 

H 12 (E) = -e^l\C l2 {E)-iu X2 {E))e- mE)/2 
1 B 2 1 

= ^- + -mQ 2 E 2 - ie, (B7) 

2m BE 2 2 ' v ; 

where m = 1/(2De) = 1/(2QA 2 ). Two complex energy variables are defined by E = E + 

i f^y-) and e = e + i — Since H± 2 is a simple harmonic oscillator Hamiltonian with 

respect to the complex energy variable E, we can calculate its propagator, (E'le - ^ 12 *!^)-^ 

Then, we can express Gi 2 (E,t\E ) in Eq. ( |B2"D explicitly, 



G l2 {E,t\E ) = =exp 
/4vrA 2 sinh(fit) 



— — — — — (e nt E 2 - e~ m £ 2 - 2££ ) 
4A 2 sinh(fit) V u / 

(B8) 



/2A 2 ,(^ + ^o)\ , , 
' - i- — — — tanh \ — )+ m 



\n 2 ft J V 2 

Straightforward Gaussian integrations over the space coordinates in Eq. ( |H6|) result in 
Eq. (pp. 

Depending on the timescale of the bath, Q, and the amount of the polarization energy 



fluctuation, A = y/2\ksT, the nonadiabatic reaction rates in Eq. ( |3.8p have two different 
limits.^! In the static bath limit (Q <C A), we take the short time limit Clt <C 1 of g(t) 
in Eq. ( |3.9|) , g(t) ~ |A 2 t 2 , then the nonadiabatic transition rate yields the Marcus rate 
given in Eq. ( |3.10|) after a Gaussian integration. The maxima of k^ A and k^ A are separated 
by 2A in this case and this is analogous to the Stokes shift observed in condensed phase 
spectroscopy.BH When the dynamics of bath degrees of freedom is very fast (Q ^> A), we 
can take the long time limit fit 3> 1 of g(t), g(t) ~ (A 2 /fl — i\)t, which leads to, 



^NA ~ ^NA ~ 2 i 2 ' (^9) 

e z + 7 Z 



where 7 = A 2 /fl. The Stokes shift between k^ A and k^ A has now disappeared as a result 
of fast fluctuation of the solvent polarization energy. 

The solvent diffusion rate is given in terms of the population Green's function Gu(E, t\Eo) 
at the crossing point E = E Q = e in Eq. ( |3.12|) . Since Ca is given by the Fokker-Planck 
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operator with a displaced harmonic potential for each electronic state, the population Green's 
functions Gu are easily obtainedllil 



Gu{E, t\E ) 



2ttA 2 (1 



-2^ 



cxp 



(E±X-(E ± A)e 
2A 2 (1 - e- 2m ) 



-nt\2 



(BIO) 



with +(— ) sign for i = 1(2), and the solvent diffusion rates can be calculated. 

In the cases of high and low barrier limits, the solvent diffusion rate given in an integral 
form in Eq. (|3.12|) can be approximately evaluated, 

. 2 



Q 
In 2 



2 
In 2 



3 In 2 



i\2 



+ ■ ■ 



pE l h < 1 
f3E l h > 1 



(Bll) 



where is the forward (i = 1) or backward (i = 2) reduced barrier height, E^ 
( e ±A) 2 /(4A). 



The Zusman's solution given in Ref. [Tl] amounts to two approximations made in the 
nonadiabatic solution. First, it is assumed that electron transfer occurs only when the 
solvation polarization energy E matches with the energy bias e, which yields, 



K(E) ss 2ttV 2 5(E-€). 



(B12) 



In this case, Eq. ( |3.2|) is further reduced to the localized transition model, and the nonadi- 
abatic transition rates k^ A reduces to the original Marcus expressions given in Eq. ( [3.10D . 
The second approximation made by Zusman is that the solvent diffusion rates are approx- 
imated by their high-barrier limits given in Eq. (|B11| ) for (3 El ^> 1, which can be written 

as 



(B13) 
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FIGURES 

FIG. 1. Eigenvalue diagrams are shown for three different values of the energy bias : (a) 
P\e\ =0 (normal regime), (b) f3\e\ = 10 (activationless regime), and (c) j3\e\ = 20 (inverted regime) 
as the coupling constant is increased. Parameters are (3\ = 10 and f3Q, = 1. 

FIG. 2. Real and imaginary parts of the eigenvalues are shown for the same three different 
values as in Fig. : (a) = 0, (b) f3\e\ = 10, and (c) 0\e\ = 20. Parameters are the same as 
those used in Fig. |l[ The real parts of the eigenvalues show bifurcation and coalescence behaviors 
as the coupling constant increases, and the imaginary parts of them asymptotically follow the Rabi 
oscillation frequency in the strong coupling limit. 

FIG. 3. Comparison of the eigenvalue solution (filled circle) with other theoretical results; the 
nonadiabatic rate (Eq. (jOf), solid line), the nonadiabatic-diffusion rate (Eq. (^), long dashed 
line), the adiabatic barrier crossing rates calculated from the adiabatic eigenvalue (Eq. ( |3.18 ), 
dot-dash line) and from the stationary- flux (Eq. (|3.19[ ), short dashed line) Solvent relaxation rates 
are chosen as (a) (30, = 0.05 (slow bath) and (b) = 1 (fast bath), respectively. Other param- 
eters are (3\ = 10 and j3e = 0. The nonadiabatic eigenvalue solution agrees with the nonadia- 
batic-diffusion rate in the weak coupling case, and the agreement with the adiabatic rates at the 
strong coupling case is qualitative. 

FIG. 4. Transition between the incoherent and coherent regimes in the symmetric electron 
transfer case is manifested in the calculation of the lowest eigenvalues. When the coupling constant 
is about half of the reorganization energy {(3\ = 10), the higher excited states cross with the first 
excited state. Symbols lying on top of the each other represent the same real parts of the complex 
conjugate eigenvalues. Other parameters are 0Q, = 1 and (3\e\ =0. 

FIG. 5. A phase diagram for the transition between the incoherent and coherent regimes in 
the asymmetric reaction case. Values of V c where two lowest nonzero eigenvalues cross (dot-dash 
line) are compared with an analytical estimate in Eq. (^^) (dashed line). Parameters are chosen 
Pfl = 0.01 and /9A = 10. 
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FIG. 6. Eigenvalue solution(filled circle) as well as the Marcus rate, Eq. ( 3.11| ), and the 



nonadiabatic transition rate, Eq. (|3.8| ), weak coupling regime are plotted a a function of the energy 
bias in the slow bath ((a), /3Q = 0.1) and the fast bath ((b), (3CI = 1) cases. Other parameters are 
j3V = 0.01. The dynamical bath effect is evident in the fast bath case in Fig. ||(b) 

FIG. 7. Dynamical bath effect is estimated by the amount of the peak shift of the forward 
nonadiabatic transition rate versus the energy bias curve from its value in the static bath case 
(f2 = 0) for various values of solvent relaxation rates, (a) When the bath is static, the nonadiabatic 
transition rate shows a maximum at |e| = A. As the bath relaxation increases, the maximum 
position of the nonadiabatic transition rate shifts into a smaller value of the energy bias, (b) The 
peak shift is plotted as a function of the solvent relaxation rate. 

FIG. 8. Time evolution of populations in the donor and acceptor states are shown for several 
values of the coupling constant. The initial wavepacket is chosen as the equilibrium distribution 
in the donor state. In weak coupling cases, (a) and (b), the population relaxes to the equilibrium 
exponentially, while in strong coupling cases, (c) and (d), the electronic coherence between the 
donor and the acceptor states is observed. 

FIG. 9. Breakdown of the positivity in the nonadiabatic diffusion equation in the case of a 
large reorganization energy (/3X 3> 1). A nonzero, real eigenvalue becomes negative in the case of 
a large reorganization energy (/3A > 12). Other parameters are j3V = 0.1 and f3Q = 1. 
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Fig. 1(a) Jung & Cao 
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Fig. 1(b) Jung & Cao 
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Fig. 1(c) Jung & Cao 
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Fig. 2 Jung & Cao 




Fig. 3(a) Jung & Cao 
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Fig. 3(b) Jung & Cao 
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Fig. 4 Jung & Cao 
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Fig. 5 Jung & Cao 
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Fig. 6 Jung & Cao 
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Fig. 7 Jung & Cao 
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Fig. 8 Jung & Cao 
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Fig. 9 Jung & Cao 
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